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Coherent structures in Dissipative Particle Dynamics simulations of the transition to 

turbulence in compressible shear flows 

Jan-Willem van de Meent, Alexander Morozov,* EUak Somfai,^ Eric Sultan, and Wim van Saarloos 
Instituut-Lorentz, University of Leiden, Postbus 9506, 2300 RA Leiden, The Netherlands 

We present simulations of coherent structures in compressible flows near the transition to turbu- 
lence using the Dissipative Particle Dynamics (DPD) method. The structures we find are remarkably 
consistent with experimental observations and DNS simulations of incompressible flows, despite a 
difference in Mach number of several orders of magnitude. The bifurcation from the laminar flow 
is bistable and shifts to higher Reynolds numbers when the fluid becomes more compressible. This 
work underlines the robustness of coherent structures in the transition to turbulence and illustrates 
the ability of particle-based methods to reproduce complex non-linear instabilities. 
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The transition to turbulence in parallel shear flows like 
pressure-driven channel flow or flow in a pipe is one of the 
classic problems of fluid mechanics. Until recently even 
predicting the correct order of magnitude for the transi- 
tional Reynolds number Re was problematic. This situ- 
ation has changed with the discovery of exact nonlinear 
solutions of the Navier-Stokes equations [1-6] . These so- 
lutions are dominated by streaks and streamwise vortices 
- low-dimensional coherent structures observed experi- 
mentally in wall-shear flows [7] which are generated via 
the self-sustaining process (SSP) proposed by Waleffe [8]. 
In the SSP, the counter-rotating quasi-streamwise vor- 
tices redistribute momentum along the wall-normal axis, 
creating spanwise modulations of the streamwise veloc- 
ity, known as streaks. The streaks in turn are subject to a 
Kelvin-Helmholtz-like instability due to the large veloc- 
ity gradient across their surface. Nonlinear interactions 
between the instability modes couple back to the original 
streamwise vortices thus closing the cycle. Even though 
the exact solutions are linearly unstable, they have been 
shown to control the transition to turbulence and turbu- 
lent dynamics at moderate Re [6, 9-14]. This scenario 
[15, 16] has emerged from a combination of intricate ex- 
periments [9, 10], large-scale numerical studies [1-6, 11- 
14] and model equations studies [8, 17]. 

In this Letter we study the robustness of the coher- 
ent structures and the self-sustaining process at the on- 
set of turbulence in compressible flows using a particle 
based method, the so-called Dissipative Particle Dynam- 
ics or DPD simulation method [18, 19]. Such a simu- 
lation method represents a fluid by discrete interacting 
particles whose motion converges to hydrodynamic be- 
havior on length-scales larger than the typical interpar- 
ticle distance. Our results not only illustrate the surpris- 
ing robustness of the coherent structures [7, 15, 16] to 
thermal fluctuations and compressibility effects, but also 
show that studies of the transition to turbulence provide 
a very attractive and informative testbed for assessing 
the strengths and weaknesses of particle based simulation 
methods. Lattice-Boltzmann methods [20] have already 
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FIG. 1: A pair of streaks in compressible plane Couette flow at 
Re = 1400 and Ma = 8. Color contours denote the streamwise 
deviation from the laminar flow, vector fields denote the av- 
erage in-plane motion, (a) Streamwise-direction view, show- 
ing s-averaged velocities, (b) Flow gradient-direction view at 
y = Q. The fast streak shows a clear sinusoidal inflection. 



been successfully applied in supercomputer turbulence 
studies, like flow past a car [21]; we demonstrate here that 
the transition to turbulence can be studied effectively on 
a regular single node computer and that the DPD re- 
sults throw new light on the robustness of the transition 
mechanism as well as on the simulation method itself. 

DPD is a well-tested and documented [18, 19] off- 
lattice method to simulate the Navier-Stokes equations. 
Its popularity is partly due to the ease of extending it 
to multiphase and viscoelastic flows. A limitation of the 
method which is not often stressed but which will come 
to the foreground here is that particle interactions have 
intrinsic timescales such that a DPD fluid is highly com- 
pressible. 

For plane Couette and pipe flows we find that at large 
enough Re, there is a hysteretic transition to a weakly- 
turbulent state dominated by coherent structures much 
like those present in DNS and experiments — see Figs.l 
and 2. As the compressibility increases, the transition 
to turbulence shifts to higher Reynolds numbers and be- 
comes less abrupt and possibly even continuous, but the 
overall features in our fluid with high compressibility and 
strong thermal fluctuations are similar to those of incom- 
pressible fluids without thermal fluctuations. 

DPD simulation method — In DPD one integrates 
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FIG. 2: Snapshots of streaky pipe flow averaged over the pipe length. The simulation box has pipe diameter d = 60 and 
periodic length I = 126. DPD parameters are (a, 7, kT) = (5, 1, 0.2), density p = 4, iV = 1425026, yielding u = 0.2, c ~ 4.2 and 
Re/Ma — 1260. (a) a 4-streak configuration at Re = 1700. (b)-(d) three difi'erent snapshots at Re = 2700. These states and 
the stochastic hopping between them are similar to those observed in experiments and in DNS of incompressible flows [9, 12]. 



Newton's equations for a system of unit-mass particles 
that represent parcels of fluid. Forces between particles 
are chosen so as to optimize the convergence to hydro- 
dynamic behavior on length scales of a few particles — a 
rationale that is similar to that of the Lattice-Boltzmann 
method [20], where microscopic fidelity is also sacrificed 
to obtain a computationally efficient representation of 
hydrodynamics beyond the lattice scale. The DPD in- 
terparticle forces are pair-wise and consist of three con- 
tributions: a soft-repulsion conservative component fc — 
a{l — r), a dissipativc component /(J = — 7(1 — r)^"((5v-f) 
that tends to reduce the difference in particle velocities 
5v, and a stochastic component /^ = ct (1 — r^^'^xit)- 
The constants a, 7 and a define the amplitude of each 
of the components; r is the distance between the parti- 
cles and f is the unit vector in the direction of f. The 
range of interaction is customarily set to 1 . A Gaussian- 
distributed random variable x with unit variance de- 
fines the evolution of random interactions. The am- 
plitude a and the form of the dissipativc and random 
forces are chosen so that a DPD fluid at rest satisfies 
the fluctuation-dissipation theorem with temperature T: 
0-2 = {2'jkT/At). The absorption of a factor y/l/At into 
a is necessary to converge properly to a continuum limit 
for small timesteps [19]. An important though seldom 
stressed point of DPD is that in order to converge quickly 
to hydrodynamic length-scales upon coarse graining, the 
three force components have to be of roughly the same 
size. This effectively limits the parameters {a,j,kT) to 
the range 0.1-10 in DPD units, and implies that density 
fluctuations, viscous interaction and thermal diffusion all 
take place on similar timescales. On the particle scale, 
DPD thus models a compressible and hot sluggish fluid. 
In our simulations, we focus on the effect of the compress- 
ibility on the transition to turbulence. One should keep 
in mind, however, that the thermal fluctuations can also 
be important — even though the thermal velocities Wth 
are much smaller than the flow velocity scale U (in our 
case typically Vth/U ^ 10~^), they determine a natural 



Reynolds number above which the flow will become un- 
stable even without external perturbations, because the 
thermal fluctuations arc sufficient to destabilize the flow 
due to the subcritical nature of the instability. 

Flow domains — Simulations were carried out in two 
classical geometries: flow between two plates y = zth slid- 
ing with opposite velocities ±t/ along each other (plane 
Couette flow), and pressure-driven flow in a pipe created 
by applying a constant force in the flow direction to ev- 
ery particle. The simulation box for the plane Couette 
geometry has dimensions 60 x 20 x 40 in the streamwise 
(velocity) x, gradient direction y, and spanwise direction 
z, and is periodic in the x and z directions. The DPD pa- 
rameters are {a,j,kT) = (8,1,0.35) unless stated other- 
wise. The density is p = 4, corresponding to A^ = 192000 
particles in the system. The parameter values for pipe 
flow are similar and given in the caption of Fig. 2. 

To impose the no-slip boundary conditions at the walls, 
we employ the method introduced by Revcnga et al. 
[22, 23] for 2D simulations. Here, the walls arc modelled 
by an immobile continuum DPD medium of uniform den- 
sity which interacts with the bulk particles (see [24] for 
details) . We have checked that in all our runs the velocity 
difference between a wall and the flrst layer of particles 
close to the wall is indeed negligibly small. 

Reynolds and Mach numbers — In our simulations, we 
characterize the flow by two parameters: the Reynolds 
number, deflned as Re — hU /v for Couette flow (for pipe 
flow Re = dU /v, where U is the averaged streamwise ve- 
locity and d the pipe diameter), and the Mach number 
Ma ~ U /c, where c is the sound velocity. We empiri- 
cally obtain c by measuring the speed of propagation of 
a density pulse directly, yielding c = 4.4 for the param- 
eter set mentioned above. The sound velocity can also 
be estimated by c ~ \Jd-pjdp\T based on the virial ex- 
pansion for the pressure p{p,T), which typically yields 
estimates within 20% of our direct measurements. The 
kinematic viscosity was measured from the laminar shear 
start-up and was found to be i^ = 0.23. Since we vary 




FIG. 3: Mean streamwise velocity U{y) and density p{y) for 
the flow snapshot shown in Fig. 1 (solid), with the initial 
t = state (dotted) shown for comparison, (a) The mean 
velocity U{y) develops a sinusoidal deviation from the linear 
profile as coherent structures develop, (b) As a result of the 
compressibility of the fluid, the density proflle p{y) shows a 
slight bulge at the center of the cell. 



U keeping other parameters fixed, we report our data in 
terms of the ratio Re/Ma = hc/i', which is of order 90- 
220 (Ma ~ 2 — 5). Thus, we are sampling a different 
parameter range than accessible in experiments [9, 10] 
and theory [1-6, 11-14]. 

Results for plane Couette flow — To quantify the in- 
fluence of the compressibility on the transition to tur- 
bulence, we focus on the plane Couette geometry. The 
qualitative features of the velocity field arc in agreement 
with the SSP predictions: Fig. la shows the stream- 
wise vortices and streaks at Re = 1400. Moreover, as 
Fig. lb shows, the streaks have a slight sinusoidal mod- 
ulation in the streamwise-spanwise plane, in agreement 
with the SSP scenario that a Kelvin-Hclmholtz-typc in- 
stability transfers energy back into the vortices. Finally, 
the mean velocity profile of Fig. 3a strongly deviates from 
the laminar profile in the way typical for turbulent flow. 
The mean density of the system, plotted in Fig. 3b, devi- 
ates slightly from its equilibrium value, which emphasizes 
the compressible nature of our fluid. 

The streaky proflle in plane Couette flow tends to 
have well-defined modes exhibiting an integer number 
of streak and vortex pairs that are fairly persistent in 
time. The bifurcation point and the dominant threshold 
mode depend on the Mach number. For low Re/Ma, the 
first mode appearing is typically a one-pair streak-vortex 
configuration as shown in Fig. 1, while for higher Re/Ma 
our simulations exhibit a bifurcation towards a two-pair 
configuration. Higher order configurations appear as the 
Reynolds number is increased further and switching be- 
tween configurations becomes more frequent. 

Reproducibility of the dynamics observed close to the 
threshold allows examination of the transition in terms of 
the bifurcation diagram of Fig. 4. We plot the deviation 
A of the normalised profile U{y)/U from the laminar flow 
for a series of states initialized at regular intervals in Re, 
marked with circles. The lower curve denotes the 2-pair 
modes which bifurcate spontaneously from the laminar 
state. The upper curve denotes 1-pair modes that are 
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FIG. 4: Bistability in the bifurcation from laminar flow at 
Re/Ma = 190. Shown is the maximum deviation A of the 
mean profile U{y) from linearity, A = maxy \{U{y)/U—y/h)\. 
Approaching from the laminar state, a smooth (forward) tran- 
sition towards a 2-pair mode is observed, which becomes un- 
stable at Re > 1200. A stable 1-pair mode is found by initial 
forcing of a pair of rolls. Decreasing the driving rate results 
in a jump (subcritical) transition back to the laminar state. 



created by initially driving all the particles with an ad- 
ditional external force term with the desired symmetry. 
The vortices that arc thus created either persist or relam- 
inarize after the forcing is turned off. The open circles 
on this branch denote states obtained this way. On both 
branches, the states can be traced smoothly by adiabati- 
cally increasing/decreasing the driving rate, thereby pro- 
ducing the continuous curves. The fact that the segments 
connect perfectly shows that the amplitude of turbulent 
perturbations is a well-defined function of Re and that 
the adjustment of U can be considered adiabatic. 

Contrary to the results for incompressible flows, the 
bifurcation observed in the lower curve of Fig. 4 appears 
continuous in Re rather than jump-like. At this point 
we are not sure which property underlies this qualitative 
difference in the onset dynamics — compressibility, flnite 
temperature T or flnite size. One possibility is that the 
2-pair mode has an instability threshold smaller than the 
typical fluctuating velocities while the resulting nonlinear 
branch lies so low that the transition appears continuous 
in Re. At the same time, the forced 1-pair state has 
presumably a higher instability threshold (since its up- 
per branch is higher than the upper branch of the 2-pair 
mode) and when tracked down along the upper branch, 
still exhibits a jump back to the laminar flow. 

We now proceed to examine how these features change 
when we change the compressibility and hence the Mach 
number Ma by tuning the repulsive force strength a. 
Fig. 5 shows 7 upper branches for a range of Ma. The 
initial states are created by seeding 1-pair modes at reg- 
ularly spaced intervals in Re, followed by a slow adjust- 
ment of the conservative force strength a towards the 
values a = 3,4,5,6,8,10,13. The corresponding sound 
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FIG. 5: Relaminarization in a series of runs with increasing 
compressibility. The box size is 48 x 16 x 32. One-pair modes 
are initialized at regular intervals in Re, with the sound veloc- 
ities ranging from c ~ 5.5 (red) to c ~ 2.8 (blue) and tracked 
down by slowly decreasing the driving rate. The resulting 
curves show that the compressibility decreases the amplitude 
of the turbulent state and the abruptness of the transition. 



velocities are c ~ 2.8, 3.2, 3.5, 3.8, 4.4, 4.9, 5.5, and so for 
fixed driving velocity the Mach numbers of these runs 
varies by a factor of 2. Subsequent interpolation of the 
states by adiabatic adjustment of the driving rate then 
results in the curves shown in the figure. Clearly there 
is an unmistakable suppression of the turbulence as the 
compressibility increases: the turbulent amplitude de- 
creases and the onset shifts to larger Re. In addition the 
bistable jump behavior changes gradually into an appar- 
ently smooth transition as the compressibility increases. 

Results for pipe flow — Fig. 2 shows snapshots of our 
results for pipe flow for two Reynolds numbers Re = 1700 
and Re = 2700. These results are consistent with the SSP 
scenario and are similar to those found in experiments 
and simulations of incompressible flows [9, 10, 12, 15, 16]. 
They also resemble exact solutions of the incompressible 
Navicr-Stokes equations [4, 5]. Streaks, visible as colored 
contours where the downstream velocity is higher (red) or 
lower (blue) than average, are stabilized by the stream- 
wise vortices (vectors). As in incompressible flows [12], 
we observe dynamical transitions between these states. 
We leave a detailed study of this process for the future. 

Conclusion — DPD simulations reproduce the qualita- 
tive features of the exact coherent structures in remark- 
able detail. Considering the large degree of compress- 
ibility in DPD fluids, the fact that the phenomenology 
differs so little constitutes new support for the SSP [8] 
as the scenario for organization of turbulence at moder- 
ate Reynolds numbers. The dependence of turbulence 
amplitudes on Ma in the DPD fluid provides clear evi- 
dence for a suppressing effect of the compressibility on 
the transition to turbulence, with apparently a crossover 
to a continuous transition in a fluctuating fluid. 

Further potential of the method lies in its flexibility in 
incorporating interactions between fluid elements. DPD 



is very easy to program, and our simulations with over 
10^ particles are quite feasible on a single node computer. 
Turbulence in multiphase and viscoelastic flows [25] , and 
in other complex fluids, clearly seems within reach. 

We thank Bruno Eckhardt and Pep Espanol for valu- 
able discussions, the EU network PHYNECS and Dutch 
science foundations NWO and FOM for support, and the 
national computer center SARA for computer time. 
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